library(MASS)
library(Zelig)

ICCdata <- read.csv(file="ICC_Zelig.csv",head=TRUE,sep=",")
attach(ICCdata)
names(ICCdata)


# Coxph SD Model 1
z.out.SDm1 <- zelig(Surv(start_q,stop_q,ratification) ~ cw_90+dem_med_hi+dem_med_hi_cw90+logmilper1+logpk+corrected_regrat100+humrightstreat+extrconflict_ongoing+const__amend__required_for_ratif+icc_elected_officials+ldrs+brit__leg__herit_1,
	model = "coxph", data = ICCdata)
z.out.SDm1

# Coxph CC Model 3 (with continuous BDeaths)
z.out.CCm3 <- zelig(Surv(start_q,stop_q,ratification) ~ dead9097+dem_med_hi+dem_dead9097+logmilper1+logpk+corrected_regrat100+humrightstreat+extrconflict_ongoing+const__amend__required_for_ratif+icc_elected_officials+ldrs+brit__leg__herit_1,
	model = "coxph", data = ICCdata)
z.out.CCm3


se3dem<-ICCdata[0:6,]
se3dem$dead9097<-as.matrix(seq(0,10000,by=2000))
se3dem$dem_dead9097<-as.matrix(seq(0,10000,by=2000))
se3dem$dem_med_hi<-as.matrix(rep(1,6), nrow=6, ncol=1)
se3dem$logmilper1<-as.matrix(rep(mean(ICCdata$logmilper1,na.rm=TRUE),6), nrow=6, ncol=1)
se3dem$logpk<-as.matrix(rep(mean(ICCdata$logpk,na.rm=TRUE),6), nrow=6, ncol=1)
se3dem$corrected_regrat100<-as.matrix(rep(mean(ICCdata$corrected_regrat100,na.rm=TRUE),6), nrow=6, ncol=1)
se3dem$humrightstreat<-as.matrix(rep(mean(ICCdata$humrightstreat,na.rm=TRUE),6), nrow=6, ncol=1)
se3dem$extrconflict_ongoing<-as.matrix(rep(mean(ICCdata$extrconflict_ongoing,na.rm=TRUE),6), nrow=6, ncol=1)
se3dem$const__amend__required_for_ratif<-as.matrix(rep(mean(ICCdata$const__amend__required_for_ratif,na.rm=TRUE),6), nrow=6, ncol=1)
se3dem$icc_elected_officials<-as.matrix(rep(mean(ICCdata$icc_elected_officials,na.rm=TRUE),6), nrow=6, ncol=1)
se3dem$ldrs<-as.matrix(rep(mean(ICCdata$ldrs,na.rm=TRUE),6), nrow=6, ncol=1)
se3dem$brit__leg__herit_1<-as.matrix(rep(mean(ICCdata$brit__leg__herit_1,na.rm=TRUE),6), nrow=6, ncol=1)

se3ndem<-se3dem
se3ndem$dem_dead9097<-as.matrix(rep(0,6), nrow=6, ncol=1)
se3ndem$dem_med_hi<-as.matrix(rep(0,6), nrow=6, ncol=1)

se3dembd0<-as.data.frame(se3dem[1,])
se3ndembd0<-as.data.frame(se3ndem[1,])
se3dembd2k<-as.data.frame(se3dem[4,])
se3ndembd2k<-as.data.frame(se3ndem[4,])

par(mfrow=c(2,2),oma=c(2,1,2,1))
plot(survfit(z.out.CCm3, newdata=se3dembd0), ylab="survival proportion", xlab="quarter", main="Democracy, 0 Deaths")
plot(survfit(z.out.CCm3, newdata=se3ndembd0), ylab="survival proportion", xlab="quarter", main="Non-democracy, 0 Deaths")
plot(survfit(z.out.CCm3, newdata=se3dembd2k), ylab="survival proportion", xlab="quarter", main="Democracy, 2,000 Deaths")
plot(survfit(z.out.CCm3, newdata=se3ndembd2k), ylab="survival proportion", xlab="quarter", main="Non-democracy, 2,000 Deaths")
title("Fig 6. Survival Curves for Ratification, by Regime Type and Conflict",
	 outer=TRUE, cex=.6)
mtext("Survival curves holding other variables at sample means, using results from Model 3.",
	 SOUTH<-1, line=0.1, adj=0, cex=.6, outer=TRUE)
mtext("Dotted lines are 95% confidence intervals.",
	 SOUTH<-1, line=0.6, adj=0, cex=.6, outer=TRUE)
